subroutine vmd_init(u0, v0, xvec, zvec, vvec, kXvec, kZvec, pvec)

  use dimensions, only : jx, jz, jv, jDirection, &   ! domain size in r-space, parameters
                         nkx, nkz, np,           &   ! limits      in k-space, set below
                         dx,  dz,  dv,           &   ! resolution  in r-space, set below
                         dkx, dkz, dp                ! resolution  in k-space, set below

  use parameters, only : phi0, kz, vPhase,  Dxxxx, Dvvvv, &            ! real,     phys
                         ny,  iWaveFrame, iFixPump, iModx, iModz, &    ! integer,  phys
                         ModAmp,  Lx, SpeckleWidth, u, &               ! real,     phys
                         restart, nt, &                                ! integer,  num
                         Vmax, dt, KxFilter, &                         ! real,     num 
                         iself, iskip, idfMod, &                       ! integer,  num 
                         iE0Modulation, iBcPerpPeriodic, iFourier, &   ! integer,  num
                         pi, iunit, isign                              ! constants


  implicit none

  real, intent(out)    :: u0, v0
  real, intent(out)    :: xvec(jx),  zvec(jz),  vvec(jv)
  real, intent(out)    :: kXvec(jx), kZvec(jz), pvec(jv)


!-- local variables ---

  integer              :: ix, iz, iv       ! counters in r-space
  integer              :: ikx, ikz, ip     ! counters in k-space
  real                 :: phi00
  real                 :: Omega0, DeltaOmega, dOmegadPhi



      Omega0(kz) = 1.106 + 1.969*(kz-0.25) +3.16*(kz-0.25)**2 - 13.0*(kz-0.25)**3  & 
             + 12.11*(kz-0.25)**4 + 1.068*(   (0.534-kz)**0.5  -  (0.534-0.25)**0.5  )   ! accurate for 0.25<kz<0.53
 
      DeltaOmega(kz,phi0) = -(4.366 - 20.86*kz + 28.09* kz**2) * ( exp(-0.5 / kz**2) / kz**2 ) * phi0**0.5  & 
                    - 343 * kz**4 * exp( -0.5 / kz**2 ) * phi0**1.5
                    !  for   0.275 < kz < 0.425
      dOmegadPhi(kz,phi0) = -(4.366 - 20.86*kz + 28.09* kz**2) * ( exp(-0.5 / kz**2) / kz**2 ) * 0.5 / phi0**0.5  & 
                    - 343 * kz**4 * exp( -0.5 / kz**2 ) * 1.5 * phi0**0.5
                    !  for   0.275 < kz < 0.425

!--  reading file with physical parameters -------------------------------------------

      open(1,file="inphysical")
      read(1,*) phi0, kz, vPhase, ny, Dxxxx, Dvvvv, iWaveFrame, iFixPump, iModx, iModz, ModAmp, Lx, SpeckleWidth, u
      close(1) ! ny is the system length in units of 2pi/kz

!--  setting spectral resolution ------------- ---------------------------------------

      dkz = kz/ny
      dkx = 2*pi / Lx
      
      v0 = -vPhase
      u0 = u


!--  echo parameters into output file  -----------------------------------------------


      open(5,file='output')

      write(5,*)'cutoff uperp in x'
      write(5,*)"code Vlasov 2D 39.f95"; write(5,*)
      write(5,*)"uPerp has periodic BC"; write(5,*)
      write(5,*)"jx=", jx, "   jz=",jz,"   jv=",jv; write(5,*)
      write(5,*)"Filtered as:  y = (xvec(ix)-Lx/4) / (Lx/03)"; write(5,*)

            
      if( iFixPump==1)then
        write(5,*)'Electric field at kz is fixed'
      else
        write(5,*)'Electric field at kz is NOT fixed'
      endif

      write(5,*)
      write(5,*)'phi0=',phi0,'    kz=',kz,' vPhase=',vPhase,'   Lx=', Lx,'     SpeckleWidth=',SpeckleWidth   
      write(5,*)'u=', u ;
     
      if( iWaveFrame==1) then
         write(5,*)'in wave frame'
      else
         write(5,*)'in lab frame'
      endif
      
!      dkz=kz/ny;  dkx = 2*pi / Lx
      
      write(5,*)
      write(5,*)'Lx=', Lx,'  speckle width=',SpeckleWidth, '  dkx/kz=', dkx/kz
      write(5,*)'actual speckle with =', SpeckleWidth, ' x 0.37233 =',0.37233*SpeckleWidth
      write(5,*)'Transit time damping rate = u/(actual speckle with) =', u / (0.37233*SpeckleWidth);
      write(5,*)
      
      
      if(ny>1) write(*,*) 'there are',ny,'  wavelengths in parallel direction '

!      v0=-vPhase

      
! unit of time is electron plasma frequency
! of length is electron Debye length
! therefore unit of speed on electron thermal speed.
      write(5,*)"external potential trapping width/ve = sqrt(e phi0/me) =",sqrt(phi0)
      write(5,*)'Hyper viscosities: Dxxxx=', Dxxxx,'   Dvvvv=', Dvvvv





!-- reading file with numerical parameters -------------------------------------------
      
      open(2,file='innumeric')
      read(2,*) restart,nt,Vmax,dt,iself, iskip, idfMod, iE0Modulation, iBcPerpPeriodic, KxFilter, iFourier
      close(2)

      !  iFourier=1:      Fourier with linearization in uz     
      !  iFourier = -1 :  upwind    
      !  iFourier = 0 :   centered difference
      !
      !  if restart.ne.0, then read from restart file, if iself=1 then self consistent potential
      !  if (idfMod=1) then insert initial density fluctuation as per terminal input
      !  if iBcPerpPeriodic = 1 then periodic in perpendicular direction(s)
      
      if( iFourier==1 )then
         write(5,*) 'advect in z with Fourier  (linearization:  uz = v0)'
      endif
      
      if( iFourier== -1 )then
         write(5,*) 'advect in z with upwind differencing'
      endif
      
      if( iFourier== 0 )then
         write(5,*) 'advect in z with centered differencing'
      endif
      
      if(   (iFourier.ne.1).and.(iFourier.ne.-1).and.(iFourier.ne.0)  ) then
             write(*,*) 'input variable iFourier not given an allowed value, stop execution after return'
             stop
      endif
      
      write(5,*)'KxFilter=', KxFilter, '   iFourier =', iFourier
      

!-- spatial resolution --

       if(iBcPerpPeriodic.NE.1)then
          dx = Lx/(jx-1)  !  if not periodic in x
          write(5,*) 'NON-periodic boundary condition'
       else
          dx = Lx/jx      !   if periodic in x        
          if(  (jx/2)*2 .NE. jx) &
               write(*,*) 'JX NOT EVEN but periodic boundary condition, make it so and recompile'
       endif
      
      dz = 2*pi/(jz*dkz)       

      write(5,*)'   dt= ', dt
      write(5,*)'vPhase * dt / dz = parallel Courant number =', vPhase * dt / dz
     
      write(5,*)
      write(5,*)'dx=',dx,'    dz=', dz, '     dz/dt=', dz/dt,'    v0 dt/dz=', v0 * dt / dz
            
      write(5,*)'phi0=',phi0,'  kz=',kz,'   V=',vPhase
      write(5,*)'  Lx = ', Lx, '  dx = Lx/jx=', dx, '  dkx =', dkx, '  jx=', jx

!-- dimension vectors in r-space --

      do ix=1,jx
         xvec(ix) = -Lx/2 + dx * (ix-1)
      enddo

      do iz=1,jz
         zvec(iz) = (iz-1)*dz
      enddo    

      do iv=1,jv
         vvec(iv)=-vPhase*iWaveFrame -Vmax + ((iv-1)/(jv-1.0))*2*Vmax
      end do
      
      dv=( vvec(jv)-vvec(1) )/(jv-1)
      !dv = 2*Vmax/(jv-1)



!-- dimension vectors in k-space --

      nkz=jz; nkx = jx; np=jv

      if( ny > nkz/2 -1) then
       write(*,*)"ny is too big,stop execution and return"
       stop
      end if


      kXvec(1)=1.0E-10      !	don't want to divide by zero later 
      do ikx=2,jx/2
         kXvec(ikx)=kXvec(ikx-1)+dkx
      end do
      kXvec(jx)=-dkx
      do ikx=jx-1,jx/2+1,-1
         kXvec(ikx)=kXvec(ikx+1)-dkx
      end do


      kZvec(1)=1.0E-10      !	don't want to divide by zero later
      do ikz=2,nkz/2
         kZvec(ikz)=kZvec(ikz-1)+dkz
      end do
      kZvec(nkz)=-dkz
      do ikz=nkz-1,nkz/2+1,-1
         kZvec(ikz)=kZvec(ikz+1)-dkz
      end do


      dp=2*pi/(jv*dv)
      pvec(1)=0.0
      do ip=2,jv/2
         pvec(ip)=pvec(ip-1)+dp
      end do
      pvec(jv)=-dp
      do ip=jv-1,jv/2+1,-1
         pvec(ip)=pvec(ip+1)-dp
      end do


      open(13,file='kXvec.txt'); 
      write(13,*) '% ik   kX' ; 
      do ikx=1,jx 
         write(13,*) ikx, kXvec(ikx)
      enddo
      close(13)


!-- dispertion relation ----


      if(iself==-1) write(5,*)'Ez2D = E0z'
      
      write(5,*)
      write(5,*) 'Domega=', -(4.366 - 20.86*kz + 28.09* kz**2) * ( exp(-0.5 / kz**2) / kz**2 ), &
                 ' x  phi^0.5 + ', - 343 * kz**4 * exp( -0.5 / kz**2 )  ,' x phi^1.5'


      if(iself.ne.1)then
        if(iself == -1) then
          write(5,*)'prescribed total potential =',phi0
        else ! iself=0
          write(5,*)'at wavenumber kz, non self consistent potential =',phi0, &
                    ' , but self-consistent harmonics allowed'
        endif
      else
       write(5,*)'self consistent, external potential =',phi0
      end if
 

      close(5)


!-- writing "Dispertion_Relation.txt"-------------------------------------------------
   
      open(1,file='Dispersion_Relation.txt'); 
      write(1,*)'% phi   Vphase     -0.5*phixdomega/dphi'
      do iz=1,40
          phi00 = iz * 0.5/40
          write(1,501) phi00, ( Omega0(kz) + DeltaOmega(kz,phi00 ) ) /  kz, & 
                       -0.5 * phi00 * dOmegadPhi(kz,phi00)
      end do
      close (1)


    501   format(  6400(e12.6,2x)  )

!----------------------------------------------------------------------------------


end subroutine vmd_init
